Show code
library(dplyr)
library(readr)
library(ggplot2)
library(tigris)
library(sf)
library(DT)Caitlin Uang
December 16, 2025
To acquire the national fast food count per state and county, this project will use the NaNDA: Eating and Drinking Places by Census Tract & ZCTA, 1990-2021 dataset and filter year 2019. The dataset at census tract level is used. NaNDA is the National Neighborhood Data Archive, providing counts and density of food establishments (restaurants, bars, fast food) from 1990-2021, sourced from the National Establishment Time Series (NETS) database.
# Cache shapefiles locally (faster repeat runs)
options(tigris_use_cache = TRUE)
# Load dataset and use only 2019 fast food count
eatdrink_tract <- read_csv("nanda_eatdrink_Tract20_1990-2021_01P.csv")
fastfood_raw <- eatdrink_tract |>
filter(year == 2019) |>
arrange(desc(count_fastfood)) |>
select(tract_fips20,totpop,count_fastfood)The U.S. Census Bureau American Community Survey (ACS) Table B01003 provides the total population for a given geographic area. The ACS B01003 table comes from 2019.
The NaNDA dataset uses tract FIPS (Federal Information Processing Standards) to identify each census tract within the U.S., while ACS B01003 uses GEOID, a full numeric code, to uniquely identify all administrative and statistical geographic areas. GEOIDs encompasses FIPS codes and will be used to join both datasets.
# Get FIPS from GEOID
pop_clean <- pop_2019 |>
mutate(tract_fips20 = substr(GEO_ID, nchar(GEO_ID) - 10, nchar(GEO_ID)))
# Join ACS and fast food dataset by FIPS code and filter out tract with 0 population
fastfood_tract <- fastfood_raw |>
left_join(pop_clean, by = "tract_fips20") |>
select(tract_fips20,GEO_ID,totpop,B01003_001E,count_fastfood) |>
rename(acs_pop = B01003_001E) |>
mutate(acs_pop = as.numeric(acs_pop)) |>
filter(acs_pop > 0)Let’s look at the joined dataset. Now that we have both dataset combined we can calculate the
#State look up
state_lookup <- states(year = 2019) |>
st_drop_geometry() |>
select(STATEFP, STUSPS, NAME) |>
rename(
State_Abbr = STUSPS,
State_Name = NAME
)
# Aggregate tract to county
den_fastfood_county <- fastfood_tract |>
mutate(county_fips20 = substr(tract_fips20, 1, 5)) |>
group_by(county_fips20) |>
summarize(
total_fastfood = sum(count_fastfood, na.rm = TRUE),
total_pop = sum(acs_pop, na.rm = TRUE),
den_fastfood = (total_fastfood / total_pop) * 1000,
.groups = "drop"
) |>
filter(total_pop > 0)
# Load shapefile and join
county_shapes <- counties(year = 2019, cb = TRUE) |>
mutate(GEOID = as.character(GEOID))
fastfood_map_county <- county_shapes |>
left_join(den_fastfood_county, by = c("GEOID" = "county_fips20")) |>
left_join(state_lookup, by = "STATEFP") |>
mutate(
den_fastfood_plot = ifelse(den_fastfood == 0, NA, den_fastfood)
)
# Create datatable
datatable(
fastfood_map_county |>
st_drop_geometry() |>
arrange(desc(den_fastfood)) |>
select(
County = NAME,
State = State_Abbr,
State_Name,
total_fastfood,
total_pop,
den_fastfood
),
caption = "Fast Food Restaurants per 1,000 People (County, 2019)",
options = list(pageLength = 10, scrollX = TRUE)
) |>
formatRound("total_fastfood", 0, mark = ",") |>
formatRound("den_fastfood", 2)fastfood_map_county |>
st_drop_geometry() |>
ggplot(aes(x = den_fastfood)) +
geom_histogram(
bins = 30,
fill = "#9badff",
color = "white",
linewidth = 0.3
) +
coord_cartesian(xlim = c(0, 1)) +
labs(
title = "Fast Food Restaurants per 1,000 People",
subtitle = "Distribution of Fast Food Density per County, 2019",
caption = "Source: NaNDA Eating and Drinking Places by Census Tract & ZCTA, 1990–2021\nU.S. Census Bureau ACS B01003 (2019)",
x = "Fast Food Density",
y = "# of Counties"
) +
theme(
axis.line = element_line(color = "black", linewidth = 0.4),
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.margin = margin(5, 5, 5, 5),
plot.caption = element_text(size = 10, face = "italic", hjust = 0),
panel.background = element_rect(fill = "white")
)ggplot(fastfood_map_county) +
geom_sf(aes(fill = den_fastfood), color = NA, linewidth = 0.1) +
scale_fill_gradientn(
colours = c("grey90", "#e792e7", "#5b7cff", "#0432FF"),
name = "Fast Food\nDensity",
limits = c(NA, NA),
na.value = "#ffffff",
trans = "sqrt"
) +
coord_sf(
xlim = c(-125, -66),
ylim = c(24, 50),
expand = FALSE
) +
labs(
title = "Fast Food Restaurants per 1,000 People",
subtitle = "County, 2019",
caption = "Source: NaNDA Eating and Drinking Places by Census Tract & ZCTA, 1990–2021\nU.S. Census Bureau ACS B01003 (2019)"
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.margin = margin(5, 5, 5, 5),
plot.caption = element_text(size = 10, face = "italic", hjust = 0)
)# Aggregate tract to state
fastfood_state <- fastfood_tract |>
mutate(state_fips20 = substr(tract_fips20, 1, 2)) |>
group_by(state_fips20) |>
summarize(
total_fastfood = sum(count_fastfood, na.rm = TRUE),
total_pop = sum(acs_pop, na.rm = TRUE),
den_fastfood = (total_fastfood / total_pop) * 1000,
.groups = "drop"
) |>
filter(total_pop > 0)
# Load state shapefile and join
state_shapes <- states(year = 2019, cb = TRUE) |>
mutate(STATEFP = as.character(STATEFP))
fastfood_map_state <- state_shapes |>
left_join(fastfood_state, by = c("STATEFP" = "state_fips20"))
# Creat datatable
datatable(
fastfood_map_state |>
st_drop_geometry() |>
arrange(desc(den_fastfood)) |>
slice_max(den_fastfood, n = 10) |>
select(State = NAME, den_fastfood,total_pop,total_fastfood),
rownames = FALSE,
caption = "Fast Food Restaurants per 1,000 People (State, 2019)",
options = list(
searching = FALSE,
scrollX = FALSE,
autoWidth = TRUE
)
) |>
formatRound("den_fastfood", 2)# Create heat map
ggplot(fastfood_map_state) +
geom_sf(aes(fill = den_fastfood), color = "white") +
scale_fill_gradient(
low = "#ffffff",
high = "#932092",
name = "Fast Food\nDensity", limits = c(NA,NA)
)+
coord_sf(
xlim = c(-125, -66),
ylim = c(24, 50),
expand = FALSE
) +
labs(
title = "Fast Food Restaurants per 1,000 People",
subtitle = "State, 2019",
caption = "Source: NaNDA Eating and Drinking Places by Census Tract & ZCTA, 1990–2021\nU.S. Census Bureau ACS B01003 (2019)"
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.margin = margin(5, 5, 5, 5),
plot.caption = element_text(size = 10, face = "italic", hjust = 0)
)# Rank and flag top 10 states
fastfood_map_state_topflag <- fastfood_map_state |>
arrange(desc(den_fastfood)) |>
mutate(
rank = row_number(),
den_top10 = ifelse(rank <= 10, den_fastfood, NA_real_)
)
ggplot(fastfood_map_state_topflag) +
geom_sf(aes(fill = den_top10),
color = "white",
linewidth = 0.3) +
scale_fill_distiller(
palette = "Blues", # Nice clean blues
direction = 1, # Light → dark
na.value = "#e0e0e0", # Light grey for non–top 10
name = "Fast-Food per\n1,000 People"
) +
labs(
title = "Top 10 States with the Most Fast-Food Restaurants per 1,000 People",
subtitle = "Non–Top 10 States Shown in Grey",
caption = "Source: NaNDA Tract20 (1990–2021) & U.S. Census ACS 2019"
) +
theme_void() +
theme(
legend.position = "right",
plot.title = element_text(size = 18, face = "bold"),
plot.subtitle = element_text(size = 12),
plot.caption = element_text(size = 9, hjust = 0),
plot.margin = margin(10, 10, 10, 10)
)